Metabarcoding Singapore
Phyloseq

1 Aim

2 Initialize

This file defines all the necessary libraries and variables

4 Map

4.1 Leaflet map

5 Environmental data

5.1 Per station

6 Phyloseq analysis

6.2 Break up into photosynthetic and non-photosynthetic

  • Opisthokonta (Metazoa, Fungi) are removed

Phyloseq All 
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 3000 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 3000 taxa by 8 taxonomic ranks ]

Phyloseq Eukaryotes 
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 668 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 668 taxa by 8 taxonomic ranks ]

Phyloseq Photosynthetic Eukaryotes 
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 271 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 271 taxa by 8 taxonomic ranks ]

Phyloseq Heterotrophic Eukaryotes 
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 397 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 397 taxa by 8 taxonomic ranks ]

Phyloseq Prokaryotes 
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 2077 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 2077 taxa by 8 taxonomic ranks ]

6.3 Normalize number of reads in each sample using median sequencing depth.

  • ! If there no cells do not transform, just set column to 0 function(x, t=total_hetero) (if(sum(x) > 0){ t * (x / sum(x))} else {x})

The median number of reads used for normalization of all is  60533
Phyloseq all
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 3000 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 3000 taxa by 8 taxonomic ranks ]

The median number of reads used for normalization of eukaryotes (auto+hetero) is  4735
Phyloseq eukaryotes (auto+hetero)
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 668 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 668 taxa by 8 taxonomic ranks ]

The median number of reads used for normalization of eukaryotes autotrophs is  3038
Phyloseq eukaryotes autotrophs
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 271 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 271 taxa by 8 taxonomic ranks ]

The median number of reads used for normalization of eukaryotes heterotrophs is  983
Phyloseq eukaryotes heterotrophs
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 397 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 397 taxa by 8 taxonomic ranks ]

The median number of reads used for normalization of prokaryotes is  54273
Phyloseq prokaryotes
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 2077 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 2077 taxa by 8 taxonomic ranks ]

6.4 Phyloseq files for abundant taxa

  • Remove taxa that are < 0.10 (euks) and <0.05 (proks) in any given sample
  • Normalize again…
Remove taxa in low abundance 

The median number of reads used for normalization of All is  24700
Phyloseq All
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 49 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 49 taxa by 8 taxonomic ranks ]

The median number of reads used for normalization of eukaryotes (auto+hetero) is  3359
Phyloseq eukaryotes (auto+hetero)
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 60 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 60 taxa by 8 taxonomic ranks ]

The median number of reads used for normalization of eukaryotes autotrophs is  2767
Phyloseq eukaryotes autotrophs
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 60 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 60 taxa by 8 taxonomic ranks ]

The median number of reads used for normalization of eukaryotes heterotrophs is  694
Phyloseq eukaryotes heterotrophs
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 83 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 83 taxa by 8 taxonomic ranks ]

The median number of reads used for normalization of prokaryotes is  24353
Phyloseq prokaryotes
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 44 taxa and 81 samples ]
sample_data() Sample Data:       [ 81 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 44 taxa by 8 taxonomic ranks ]

6.7 Treemaps at division and class levels

6.7.2 Do the individual treemaps

array_treemap_all <- list()
array_treemap_by_kingdom <- list()

for (one_strait in c("Singapore", "Johor")) {
    
    label <- str_c(one_strait, "-all")
    
    array_treemap_all[[label]] <- treemap_gg_dv(filter(long_all, strait == one_strait), 
        kingdom, supergroup, str_c("", one_strait))
    
    label <- str_c(one_strait, "-arch")
    array_treemap_by_kingdom[[label]] <- treemap_gg_dv(filter(long_prok, strait == 
        one_strait & kingdom == "Archaea"), division, class, str_c("A - Archaea - ", 
        one_strait))
    
    label <- str_c(one_strait, "-bact")
    array_treemap_by_kingdom[[label]] <- treemap_gg_dv(filter(long_prok, strait == 
        one_strait & kingdom == "Bacteria"), division, class, str_c("B - Bacteria - ", 
        one_strait))
    
    
    label <- str_c(one_strait, "-euks")
    array_treemap_by_kingdom[[label]] <- treemap_gg_dv(filter(long_euk, strait == 
        one_strait), division, class, str_c("C - All Eulkaryotes - ", one_strait))
    
    
    # label <- str_c(one_strait,'-photeuks') array_treemap[[label]] <-
    # treemap_gg_dv(filter(long_photo, strait == one_strait), division, class,
    # str_c('D - Photosynthetic Eukaryotes - ', one_strait ))
    
    
    # treemap_dv(filter(long_euk, strait == one_strait), c('division',
    # 'class'),'n_seq', str_c('All euks - ', one_strait ))
    # treemap_dv(filter(long_photo, strait == one_strait), c('division',
    # 'class'),'n_seq',str_c('Photo euks', one_strait ))
    # treemap_dv(filter(long_prok, strait == one_strait & kingdom=='Bacteria'),
    # c('division', 'class'),'n_seq',str_c('Bacteria - ', one_strait ))
    # treemap_dv(filter(long_prok, strait == one_strait & kingdom=='Archaea'),
    # c('division', 'class'),'n_seq',str_c('Archaea - ', one_strait ))
}

6.7.3 FIG - Arrange the different treemaps in a grid and save

TableGrob (1 x 2) "arrange": 2 grobs
              z     cells    name           grob
Singapore-all 1 (1-1,1-1) arrange gtable[layout]
Johor-all     2 (1-1,2-2) arrange gtable[layout]

TableGrob (3 x 2) "arrange": 6 grobs
               z     cells    name           grob
Singapore-arch 1 (1-1,1-1) arrange gtable[layout]
Singapore-bact 2 (2-2,1-1) arrange gtable[layout]
Singapore-euks 3 (3-3,1-1) arrange gtable[layout]
Johor-arch     4 (1-1,2-2) arrange gtable[layout]
Johor-bact     5 (2-2,2-2) arrange gtable[layout]
Johor-euks     6 (3-3,2-2) arrange gtable[layout]

6.8 Most abundant taxa

6.8.3 FIG Most abundant taxa - Save as Grid

6.14 Time series (abundant OTUs) - Division level

6.15 Time series (abundant OTUs) - Genus level for Chlorophyta

6.16 Main species for each division (Eukaryotes - Autrotrophs)

6.17 Heatmaps

6.18 NMDS

Sample removed because they were pulling the NMDS * PR2X16XS21 it has a single eukaryote (diatom bloom ) * RM13XS36 cause problem for bacteria * PR11XS25 cause problem for hetero euks * SBW11XS26 cause problem for hetero euks * SBW13XS37 cause problem for hetero euks * RM13XS36 cause problem for hetero euks

6.18.1 Define function

  • See comments inside functions (saved plots are different from displayed plots)
ps_do_nmds <- function(ps_list) {
    
    plot_array <- list()
    
    for (i in 1:3) {
        ps_nmds <- ps_list$ps[[i]]
        
        # Remove samples with no reads
        ps_nmds <- prune_samples(sample_sums(ps_nmds) > 0, ps_nmds)
        
        # Remove samples from Raffles
        ps_nmds <- subset_samples(ps_nmds, !(str_detect(strait, "Raffles Mari|Johor")))
        
        # Remove samples that caused problems (1= prok, 2=euk, 3=euk auto)
        if (i == 1) 
            {
                ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("RM13XS36")), 
                  ps_nmds)
            }  # Prokaryotes
        if (i %in% c(2, 3)) 
            {
                ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("PR2X16SXS21")), 
                  ps_nmds)
            }  # Eukaryotes
        if (i == 4) 
            {
                ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("PR11XS25", 
                  "SBW11XS26", "SBW13XS37")), ps_nmds)
            }  # Heterotrophs not used
        
        singa.ord <- ordinate(ps_nmds, "NMDS", "bray")
        
        # Fit environmental parameters
        env_var <- sample_variables(ps_nmds)
        env_matrix <- get_variable(ps_nmds, c("Chl", "Temperature", "Salinity", 
            "Phosphate", "Silicate", "DIN", "BAC"))
        env_fit <- vegan::envfit(singa.ord, env = env_matrix, perm = 999, na.rm = TRUE)
        env_arrows <- data.frame(env_fit$vectors$arrows * sqrt(env_fit$vectors$r)) %>% 
            rownames_to_column(var = "parameter")
        
        
        nmds_samples <- data.frame(singa.ord[["points"]], get_variable(ps_nmds, 
            c("monsoon", "strait", "location", "location_label"))) %>% rownames_to_column(var = "sample")
        
        # Factor to move the labels
        nudge_x <- max(nmds_samples$MDS1) * 0.08
        nudge_y <- max(nmds_samples$MDS2) * 0.08
        xy_max = max(c(nmds_samples$MDS1, nmds_samples$MDS2)) * 1.5
        xy_min = min(c(nmds_samples$MDS1, nmds_samples$MDS2)) * 1.5
        factor <- 3  # for vectors for euks
        if (i == 1) 
            {
                factor <- 1.5
            }  # for vectors for proks
        print(factor)
        
        p <- plot_ordination(ps_nmds, singa.ord, type = "samples", color = "monsoon", 
            shape = "strait", title = ps_list$title[[i]]) + geom_point(aes(shape = strait, 
            color = monsoon), size = 3.5) + scale_color_viridis(discrete = TRUE) + 
            scale_shape_manual(values = c(15, 16)) + geom_text(aes(label = location_label, 
            color = monsoon), nudge_x = nudge_x, nudge_y = nudge_y, check_overlap = TRUE, 
            size = 2) + theme_bw() + geom_segment(data = env_arrows, aes(x = 0, 
            xend = NMDS1 * factor, y = 0, yend = NMDS2 * factor), inherit.aes = FALSE, 
            arrow = arrow(length = unit(0.5, "cm")), colour = "black") + geom_text(data = env_arrows, 
            aes(x = NMDS1 * factor, y = NMDS2 * factor, label = parameter), 
            inherit.aes = FALSE, hjust = -0.2, vjust = -0.2, size = 3)
        print(singa.ord)
        plot_array[[i]] <- p
        print(p)
        
        # The following lines can be used if you want to avoid using the pjhyloseq
        # functions to plot the data.  Notes : - must use inherit.aes = FALSE to add
        # some extra layers - the saved plots have a different scale for the added
        # layer than the displaued plot can figure out ggplot()+ coord_fixed() +
        # xlim(xy_min, xy_max) + ylim(xy_min, xy_max) +
        # geom_point(data=nmds_samples, aes(x=MDS1, y=MDS2, shape=strait,
        # color=monsoon), size=5) + geom_text(data=nmds_samples, aes(x=MDS1, y=MDS2,
        # label=location_label, color=monsoon), nudge_x=nudge_x, nudge_y=nudge_y,
        # check_overlap = FALSE, size=2) + ggtitle(ps_list$title[[i]]) +
        
        p <- plot_ordination(ps_nmds, singa.ord, type = "taxa", color = "division", 
            title = ps_list$title[[i]]) + scale_color_viridis(discrete = TRUE) + 
            geom_point(size = 3) + theme_bw() + geom_segment(data = env_arrows, 
            aes(x = 0, xend = NMDS1 * factor, y = 0, yend = NMDS2 * factor), 
            inherit.aes = FALSE, arrow = arrow(length = unit(0.5, "cm")), colour = "black") + 
            geom_text(data = env_arrows, aes(x = NMDS1 * factor, y = NMDS2 * 
                factor, label = parameter), inherit.aes = FALSE, hjust = -0.2, 
                vjust = -0.2, size = 3)
        
        print(p)
        plot_array[[i + 3]] <- p
    }
    return(plot_array)
}

6.18.2 All OTUs

Square root transformation
Wisconsin double standardization
Run 0 stress 0.12 
Run 1 stress 0.12 
... New best solution
... Procrustes: rmse 0.032  max resid 0.14 
Run 2 stress 0.12 
Run 3 stress 0.12 
... Procrustes: rmse 0.003  max resid 0.013 
Run 4 stress 0.12 
... New best solution
... Procrustes: rmse 0.00014  max resid 0.00035 
... Similar to previous best
Run 5 stress 0.16 
Run 6 stress 0.12 
... New best solution
... Procrustes: rmse 5.4e-05  max resid 0.00019 
... Similar to previous best
Run 7 stress 0.12 
... Procrustes: rmse 0.0031  max resid 0.013 
Run 8 stress 0.18 
Run 9 stress 0.15 
Run 10 stress 0.18 
Run 11 stress 0.12 
Run 12 stress 0.12 
... Procrustes: rmse 0.0031  max resid 0.013 
Run 13 stress 0.12 
... Procrustes: rmse 0.00011  max resid 0.00041 
... Similar to previous best
Run 14 stress 0.12 
Run 15 stress 0.12 
Run 16 stress 0.17 
Run 17 stress 0.12 
... Procrustes: rmse 6.3e-05  max resid 0.00026 
... Similar to previous best
Run 18 stress 0.12 
... Procrustes: rmse 0.0029  max resid 0.012 
Run 19 stress 0.15 
Run 20 stress 0.12 
*** Solution reached
[1] 1.5

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.12 
Stress type 1, weak ties
Two convergent solutions found after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))' 

Square root transformation
Wisconsin double standardization
Run 0 stress 0.22 
Run 1 stress 0.22 
... Procrustes: rmse 0.018  max resid 0.08 
Run 2 stress 0.24 
Run 3 stress 0.24 
Run 4 stress 0.22 
Run 5 stress 0.23 
Run 6 stress 0.22 
... New best solution
... Procrustes: rmse 0.065  max resid 0.16 
Run 7 stress 0.26 
Run 8 stress 0.22 
... Procrustes: rmse 0.074  max resid 0.29 
Run 9 stress 0.22 
... New best solution
... Procrustes: rmse 0.047  max resid 0.15 
Run 10 stress 0.23 
Run 11 stress 0.22 
Run 12 stress 0.24 
Run 13 stress 0.22 
... Procrustes: rmse 0.034  max resid 0.14 
Run 14 stress 0.22 
Run 15 stress 0.24 
Run 16 stress 0.22 
Run 17 stress 0.23 
Run 18 stress 0.22 
... New best solution
... Procrustes: rmse 0.02  max resid 0.089 
Run 19 stress 0.25 
Run 20 stress 0.23 
*** No convergence -- monoMDS stopping criteria:
    20: stress ratio > sratmax
[1] 3

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.22 
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))' 

Square root transformation
Wisconsin double standardization
Run 0 stress 0.2 
Run 1 stress 0.19 
... New best solution
... Procrustes: rmse 0.093  max resid 0.32 
Run 2 stress 0.23 
Run 3 stress 0.19 
... New best solution
... Procrustes: rmse 0.00045  max resid 0.0019 
... Similar to previous best
Run 4 stress 0.18 
... New best solution
... Procrustes: rmse 0.054  max resid 0.2 
Run 5 stress 0.18 
... Procrustes: rmse 0.018  max resid 0.072 
Run 6 stress 0.22 
Run 7 stress 0.19 
Run 8 stress 0.22 
Run 9 stress 0.19 
Run 10 stress 0.21 
Run 11 stress 0.21 
Run 12 stress 0.18 
... Procrustes: rmse 0.00022  max resid 0.0011 
... Similar to previous best
Run 13 stress 0.19 
Run 14 stress 0.2 
Run 15 stress 0.18 
Run 16 stress 0.18 
... New best solution
... Procrustes: rmse 7.5e-05  max resid 0.00035 
... Similar to previous best
Run 17 stress 0.2 
Run 18 stress 0.18 
... Procrustes: rmse 0.00012  max resid 0.00049 
... Similar to previous best
Run 19 stress 0.18 
... Procrustes: rmse 0.018  max resid 0.072 
Run 20 stress 0.18 
... Procrustes: rmse 3.6e-05  max resid 1e-04 
... Similar to previous best
*** Solution reached
[1] 3

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.18 
Stress type 1, weak ties
Two convergent solutions found after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'

6.18.3 Abundant OTUs

Square root transformation
Wisconsin double standardization
Run 0 stress 0.17 
Run 1 stress 0.17 
... Procrustes: rmse 0.00024  max resid 0.001 
... Similar to previous best
Run 2 stress 0.17 
Run 3 stress 0.19 
Run 4 stress 0.19 
Run 5 stress 0.19 
Run 6 stress 0.22 
Run 7 stress 0.17 
... New best solution
... Procrustes: rmse 0.019  max resid 0.08 
Run 8 stress 0.18 
Run 9 stress 0.17 
Run 10 stress 0.2 
Run 11 stress 0.19 
Run 12 stress 0.2 
Run 13 stress 0.17 
... Procrustes: rmse 4.9e-05  max resid 0.00014 
... Similar to previous best
Run 14 stress 0.16 
... New best solution
... Procrustes: rmse 0.054  max resid 0.2 
Run 15 stress 0.18 
Run 16 stress 0.17 
Run 17 stress 0.18 
Run 18 stress 0.2 
Run 19 stress 0.18 
Run 20 stress 0.16 
... Procrustes: rmse 0.00012  max resid 0.00045 
... Similar to previous best
*** Solution reached
[1] 1.5

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.16 
Stress type 1, weak ties
Two convergent solutions found after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))' 

Square root transformation
Wisconsin double standardization
Run 0 stress 0.21 
Run 1 stress 0.25 
Run 2 stress 0.22 
Run 3 stress 0.21 
... New best solution
... Procrustes: rmse 0.00014  max resid 0.00043 
... Similar to previous best
Run 4 stress 0.24 
Run 5 stress 0.21 
... Procrustes: rmse 5.6e-05  max resid 0.00019 
... Similar to previous best
Run 6 stress 0.23 
Run 7 stress 0.25 
Run 8 stress 0.26 
Run 9 stress 0.24 
Run 10 stress 0.21 
... Procrustes: rmse 0.019  max resid 0.086 
Run 11 stress 0.21 
... New best solution
... Procrustes: rmse 0.00012  max resid 0.00041 
... Similar to previous best
Run 12 stress 0.22 
Run 13 stress 0.21 
... Procrustes: rmse 0.00052  max resid 0.0017 
... Similar to previous best
Run 14 stress 0.24 
Run 15 stress 0.28 
Run 16 stress 0.22 
Run 17 stress 0.22 
Run 18 stress 0.24 
Run 19 stress 0.21 
... Procrustes: rmse 0.00032  max resid 0.0011 
... Similar to previous best
Run 20 stress 0.21 
... Procrustes: rmse 9.6e-05  max resid 3e-04 
... Similar to previous best
*** Solution reached
[1] 3

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.21 
Stress type 1, weak ties
Two convergent solutions found after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))' 

Square root transformation
Wisconsin double standardization
Run 0 stress 0.21 
Run 1 stress 0.23 
Run 2 stress 0.22 
Run 3 stress 0.2 
... New best solution
... Procrustes: rmse 0.067  max resid 0.27 
Run 4 stress 0.24 
Run 5 stress 0.21 
Run 6 stress 0.24 
Run 7 stress 0.24 
Run 8 stress 0.23 
Run 9 stress 0.25 
Run 10 stress 0.23 
Run 11 stress 0.21 
Run 12 stress 0.22 
Run 13 stress 0.21 
Run 14 stress 0.2 
Run 15 stress 0.22 
Run 16 stress 0.21 
Run 17 stress 0.21 
Run 18 stress 0.23 
Run 19 stress 0.2 
Run 20 stress 0.23 
*** No convergence -- monoMDS stopping criteria:
     1: no. of iterations >= maxit
    19: stress ratio > sratmax
[1] 3

Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
Distance: bray 

Dimensions: 2 
Stress:     0.2 
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'

### MSDS graph

6.20 Differential expression (DESeq2)

Love, M.I., Huber, W. & Anders, S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550.

[1] '1.22.2'

Block 1 is need to avoid the following error > Error in estimateSizeFactorsForMatrix(counts(object), locfunc = locfunc, : every gene contains at least one zero, cannot compute log geometric means

See: https://github.com/joey711/phyloseq/issues/387

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 3000 taxa and 21 samples ]
sample_data() Sample Data:       [ 21 samples by 26 sample variables ]
tax_table()   Taxonomy Table:    [ 3000 taxa by 8 taxonomic ranks ]
         baseMean log2FoldChange lfcSE stat  pvalue    padj  kingdom
otu_0005     2579            4.8  1.48  3.3 0.00104 0.00501  Archaea
otu_0017      728            1.5  0.48  3.2 0.00156 0.00723 Bacteria
otu_0025       42            9.0  2.60  3.5 0.00054 0.00266 Bacteria
otu_0054      616            1.5  0.44  3.5 0.00052 0.00261 Bacteria
otu_0074      261            4.9  1.27  3.8 0.00012 0.00067 Bacteria
             supergroup            division            class
otu_0005 Thaumarchaeota     Nitrososphaeria Nitrosopumilales
otu_0017 Proteobacteria Gammaproteobacteria      SAR86 clade
otu_0025 Actinobacteria      Acidimicrobiia   Microtrichales
otu_0054 Proteobacteria Alphaproteobacteria      SAR11 clade
otu_0074 Proteobacteria Alphaproteobacteria      SAR11 clade
                      order                    family
otu_0005  Nitrosopumilaceae Candidatus Nitrosopumilus
otu_0017        SAR86 clade               SAR86 clade
otu_0025 Ilumatobacteraceae             Ilumatobacter
otu_0054            Clade I                  Clade Ib
otu_0074          Clade III                 Clade III
                             genus species
otu_0005 Candidatus Nitrosopumilus    <NA>
otu_0017               SAR86 clade    <NA>
otu_0025             Ilumatobacter    <NA>
otu_0054                  Clade Ib    <NA>
otu_0074                 Clade III    <NA>
 [ reached 'max' / getOption("max.print") -- omitted 1 rows ]

Daniel Vaulot, Adriana Lopes dos Santos

22 04 2019